# install packages
install.packages("statnet")
install.packages("xergm")
install.packages("texreg")
install.packages("xtable")
install.packages("ggplot2")

# set directories
setwd("~/Dropbox/ISQ_ms/5 Final/Replication")
SAVEDIR <- "~/Dropbox/ISQ_ms/5 Final/Replication"

# load libraries
detach(package:igraph, unload=T)
require("statnet")
require("xergm")
require("texreg")
require("xtable")
library(ggplot2)


################################################ Model1a ############################################################

### Load data

# Read in edgelists
el1970 <- read.csv("mod1_el1970.csv",stringsAsFactors=F)
el1975 <- read.csv("mod1_el1975.csv",stringsAsFactors=F)
el1980 <- read.csv("mod1_el1980.csv",stringsAsFactors=F)
el1985 <- read.csv("mod1_el1985.csv",stringsAsFactors=F)
el1990 <- read.csv("mod1_el1990.csv",stringsAsFactors=F)
el1995 <- read.csv("mod1_el1995.csv",stringsAsFactors=F)
el2000 <- read.csv("mod1_el2000.csv",stringsAsFactors=F)
el2005 <- read.csv("mod1_el2005.csv",stringsAsFactors=F)

# Read in dyadic covariates

# alliance networks
ally1970 <- read.csv("mod1_ally1970.csv",row.names=1)
ally1970 <- as.matrix(ally1970)
ally1975 <- read.csv("mod1_ally1975.csv",row.names=1)
ally1975 <- as.matrix(ally1975)
ally1980 <- read.csv("mod1_ally1980.csv",row.names=1)
ally1980 <- as.matrix(ally1980)
ally1985 <- read.csv("mod1_ally1985.csv",row.names=1)
ally1985 <- as.matrix(ally1985)
ally1990 <- read.csv("mod1_ally1990.csv",row.names=1)
ally1990 <- as.matrix(ally1990)
ally1995 <- read.csv("mod1_ally1995.csv",row.names=1)
ally1995 <- as.matrix(ally1995)
ally2000 <- read.csv("mod1_ally2000.csv",row.names=1)
ally2000 <- as.matrix(ally2000)
ally2005 <- read.csv("mod1_ally2003.csv",row.names=1)
ally2005 <- as.matrix(ally2005)

# contiguity matrix
cont1970 <- read.csv("mod1_cont1970.csv",row.names=1)
cont1970 <- as.matrix(cont1970)
cont1975 <- read.csv("mod1_cont1975.csv",row.names=1)
cont1975 <- as.matrix(cont1975)
cont1980 <- read.csv("mod1_cont1980.csv",row.names=1)
cont1980 <- as.matrix(cont1980)
cont1985 <- read.csv("mod1_cont1985.csv",row.names=1)
cont1985 <- as.matrix(cont1985)
cont1990 <- read.csv("mod1_cont1990.csv",row.names=1)
cont1990 <- as.matrix(cont1990)
cont1995 <- read.csv("mod1_cont1995.csv",row.names=1)
cont1995 <- as.matrix(cont1995)
cont2000 <- read.csv("mod1_cont2000.csv",row.names=1)
cont2000 <- as.matrix(cont2000)
cont2005 <- read.csv("mod1_cont2005.csv",row.names=1)
cont2005 <- as.matrix(cont2005)

# trade network
trade1970 <- read.csv("mod1_trade1970.csv",row.names=1)
trade1970 <- as.matrix(trade1970)
trade1975 <- read.csv("mod1_trade1975.csv",row.names=1)
trade1975 <- as.matrix(trade1975)
trade1980 <- read.csv("mod1_trade1980.csv",row.names=1)
trade1980 <- as.matrix(trade1980)
trade1985 <- read.csv("mod1_trade1985.csv",row.names=1)
trade1985 <- as.matrix(trade1985)
trade1990 <- read.csv("mod1_trade1990.csv",row.names=1)
trade1990 <- as.matrix(trade1990)
trade1995 <- read.csv("mod1_trade1995.csv",row.names=1)
trade1995 <- as.matrix(trade1995)
trade2000 <- read.csv("mod1_trade2000.csv",row.names=1)
trade2000 <- as.matrix(trade2000)
trade2005 <- read.csv("mod1_trade2005.csv",row.names=1)
trade2005 <- as.matrix(trade2005)

# Read in vertex covariates
vcov <- read.csv("mod1_vcov.csv",stringsAsFactors=F)
vcov1970 <- vcov[vcov$year==1970,]
vcov1975 <- vcov[vcov$year==1975,]
vcov1980 <- vcov[vcov$year==1980,]
vcov1985 <- vcov[vcov$year==1985,]
vcov1990 <- vcov[vcov$year==1990,]
vcov1995 <- vcov[vcov$year==1995,]
vcov2000 <- vcov[vcov$year==2000,]
vcov2005 <- vcov[vcov$year==2005,]

### Create network

# Initialize networks and assign labels

net1970 <- network.initialize(nrow(vcov1970))
network.vertex.names(net1970) <- vcov1970$stateabb

net1975 <- network.initialize(nrow(vcov1975))
network.vertex.names(net1975) <- vcov1975$stateabb

net1980 <- network.initialize(nrow(vcov1980))
network.vertex.names(net1980) <- vcov1980$stateabb

net1985 <- network.initialize(nrow(vcov1985))
network.vertex.names(net1985) <- vcov1985$stateabb

net1990 <- network.initialize(nrow(vcov1990))
network.vertex.names(net1990) <- vcov1990$stateabb

net1995 <- network.initialize(nrow(vcov1995))
network.vertex.names(net1995) <- vcov1995$stateabb

net2000 <- network.initialize(nrow(vcov2000))
network.vertex.names(net2000) <- vcov2000$stateabb

net2005 <- network.initialize(nrow(vcov2005))
network.vertex.names(net2005) <- vcov2005$stateabb


# Add in edges
net1970[as.matrix(el1970)] <- 1
net1975[as.matrix(el1975)] <- 1
net1980[as.matrix(el1980)] <- 1
net1985[as.matrix(el1985)] <- 1
net1990[as.matrix(el1990)] <- 1
net1995[as.matrix(el1995)] <- 1
net2000[as.matrix(el2000)] <- 1
net2005[as.matrix(el2005)] <- 1

# Define vertex attributes

# For 1970

set.vertex.attribute(net1970,"pop", vcov1970$log_pop)
set.vertex.attribute(net1970,"gdp", vcov1970$log_realgdp)
set.vertex.attribute(net1970,"gdppc", vcov1970$log_rgdppc)

set.vertex.attribute(net1970,"cinc", vcov1970$log_cinc)
set.vertex.attribute(net1970,"milex", vcov1970$log_milex)
set.vertex.attribute(net1970,"milper", vcov1970$log_milper)

set.vertex.attribute(net1970,"nw", vcov1970$nw)

set.vertex.attribute(net1970,"polity", vcov1970$polity)
set.vertex.attribute(net1970,"dem_dum", vcov1970$dem_dum)

set.vertex.attribute(net1970,"region", vcov1970$sub_region_n)
set.vertex.attribute(net1970,"igo_hq", vcov1970$igo_hq)

# For 1975

set.vertex.attribute(net1975,"pop", vcov1975$log_pop)
set.vertex.attribute(net1975,"gdp", vcov1975$log_realgdp)
set.vertex.attribute(net1975,"gdppc", vcov1975$log_rgdppc)

set.vertex.attribute(net1975,"cinc", vcov1975$log_cinc)
set.vertex.attribute(net1975,"milex", vcov1975$log_milex)
set.vertex.attribute(net1975,"milper", vcov1975$log_milper)

set.vertex.attribute(net1975,"nw", vcov1975$nw)

set.vertex.attribute(net1975,"polity", vcov1975$polity)
set.vertex.attribute(net1975,"dem_dum", vcov1975$dem_dum)

set.vertex.attribute(net1975,"region", vcov1975$sub_region_n)
set.vertex.attribute(net1975,"igo_hq", vcov1975$igo_hq)

# For 1980

set.vertex.attribute(net1980,"pop", vcov1980$log_pop)
set.vertex.attribute(net1980,"gdp", vcov1980$log_realgdp)
set.vertex.attribute(net1980,"gdppc", vcov1980$log_rgdppc)

set.vertex.attribute(net1980,"cinc", vcov1980$log_cinc)
set.vertex.attribute(net1980,"milex", vcov1980$log_milex)
set.vertex.attribute(net1980,"milper", vcov1980$log_milper)

set.vertex.attribute(net1980,"nw", vcov1980$nw)

set.vertex.attribute(net1980,"polity", vcov1980$polity)
set.vertex.attribute(net1980,"dem_dum", vcov1980$dem_dum)

set.vertex.attribute(net1980,"region", vcov1980$sub_region_n)
set.vertex.attribute(net1980,"igo_hq", vcov1980$igo_hq)


# For 1985

set.vertex.attribute(net1985,"pop", vcov1985$log_pop)
set.vertex.attribute(net1985,"gdp", vcov1985$log_realgdp)
set.vertex.attribute(net1985,"gdppc", vcov1985$log_rgdppc)

set.vertex.attribute(net1985,"cinc", vcov1985$log_cinc)
set.vertex.attribute(net1985,"milex", vcov1985$log_milex)
set.vertex.attribute(net1985,"milper", vcov1985$log_milper)

set.vertex.attribute(net1985,"nw", vcov1985$nw)

set.vertex.attribute(net1985,"polity", vcov1985$polity)
set.vertex.attribute(net1985,"dem_dum", vcov1985$dem_dum)

set.vertex.attribute(net1985,"region", vcov1985$sub_region_n)
set.vertex.attribute(net1985,"igo_hq", vcov1985$igo_hq)

# For 1990

set.vertex.attribute(net1990,"pop", vcov1990$log_pop)
set.vertex.attribute(net1990,"gdp", vcov1990$log_realgdp)
set.vertex.attribute(net1990,"gdppc", vcov1990$log_rgdppc)

set.vertex.attribute(net1990,"cinc", vcov1990$log_cinc)
set.vertex.attribute(net1990,"milex", vcov1990$log_milex)
set.vertex.attribute(net1990,"milper", vcov1990$log_milper)

set.vertex.attribute(net1990,"nw", vcov1990$nw)

set.vertex.attribute(net1990,"polity", vcov1990$polity)
set.vertex.attribute(net1990,"dem_dum", vcov1990$dem_dum)

set.vertex.attribute(net1990,"region", vcov1990$sub_region_n)
set.vertex.attribute(net1990,"igo_hq", vcov1990$igo_hq)

# For 1995

set.vertex.attribute(net1995,"pop", vcov1995$log_pop)
set.vertex.attribute(net1995,"gdp", vcov1995$log_realgdp)
set.vertex.attribute(net1995,"gdppc", vcov1995$log_rgdppc)

set.vertex.attribute(net1995,"cinc", vcov1995$log_cinc)
set.vertex.attribute(net1995,"milex", vcov1995$log_milex)
set.vertex.attribute(net1995,"milper", vcov1995$log_milper)

set.vertex.attribute(net1995,"nw", vcov1995$nw)

set.vertex.attribute(net1995,"polity", vcov1995$polity)
set.vertex.attribute(net1995,"dem_dum", vcov1995$dem_dum)

set.vertex.attribute(net1995,"region", vcov1995$sub_region_n)
set.vertex.attribute(net1995,"igo_hq", vcov1995$igo_hq)

# For 2000

set.vertex.attribute(net2000,"pop", vcov2000$log_pop)
set.vertex.attribute(net2000,"gdp", vcov2000$log_realgdp)
set.vertex.attribute(net2000,"gdppc", vcov2000$log_rgdppc)

set.vertex.attribute(net2000,"cinc", vcov2000$log_cinc)
set.vertex.attribute(net2000,"milex", vcov2000$log_milex)
set.vertex.attribute(net2000,"milper", vcov2000$log_milper)

set.vertex.attribute(net2000,"nw", vcov2000$nw)

set.vertex.attribute(net2000,"polity", vcov2000$polity)
set.vertex.attribute(net2000,"dem_dum", vcov2000$dem_dum)

set.vertex.attribute(net2000,"region", vcov2000$sub_region_n)
set.vertex.attribute(net2000,"igo_hq", vcov2000$igo_hq)

# For 2005

set.vertex.attribute(net2005,"pop", vcov2005$log_pop)
set.vertex.attribute(net2005,"gdp", vcov2005$log_realgdp)
set.vertex.attribute(net2005,"gdppc", vcov2005$log_rgdppc)

set.vertex.attribute(net2005,"cinc", vcov2005$log_cinc)
set.vertex.attribute(net2005,"milex", vcov2005$log_milex)
set.vertex.attribute(net2005,"milper", vcov2005$log_milper)

set.vertex.attribute(net2005,"nw", vcov2005$nw)

set.vertex.attribute(net2005,"polity", vcov2005$polity)
set.vertex.attribute(net2005,"dem_dum", vcov2005$dem_dum)

set.vertex.attribute(net2005,"region", vcov2005$sub_region_n)
set.vertex.attribute(net2005,"igo_hq", vcov2005$igo_hq)


# Check networks
net1970
net1975
net1980
net1985
net1990
net1995
net2000
net2005


### For the 1970-2005 period

# t=1 -> 1970
# t=2 -> 1975
# t=3 -> 1980
# t=4 -> 1985
# t=5 -> 1990
# t=6 -> 1995
# t=7 -> 2000
# t=8 -> 2005

# Combine dependent networks 

all_dip_nets <- list(as.matrix(net1970), as.matrix(net1975), as.matrix(net1980), as.matrix(net1985),
                     as.matrix(net1990), as.matrix(net1995), as.matrix(net2000), as.matrix(net2005))

# Combine dyadic covariates

all_ally_nets <- list(ally1970, ally1975, ally1980, ally1985, ally1990, ally1995, ally2000)
all_cont_nets <- list(cont1970, cont1975, cont1980, cont1985, cont1990, cont1995, cont2000)
all_trade_nets <- list(trade1970, trade1975, trade1980, trade1985, trade1990, trade1995, trade2000)


# check objects
class(all_dip_nets)
length(all_dip_nets)

# Use the preprocess function on the outcome network
# Tell it what missing values and structural zeros (instances in which tie impossible) look like

# estimation (dep nets) starts at t=2 and ends at t=8
all_dep <- preprocess(all_dip_nets,
                      lag = TRUE, covariate = FALSE, na = NA,
                      na.method = "fillmode", structzero = 10,
                      structzero.method = "remove") 

# Create lagged dependent network
all_mem <- preprocess(all_dip_nets, 
                      lag = TRUE, covariate = TRUE, memory = "autoregression",
                      na = NA, na.method = "fillmode", structzero = 10,
                      structzero.method = "remove")

# Check that everything looks right
length(all_dep)
sapply(all_dep, dim)
sapply(all_mem, dim)
# dim of the mem term need to be the same (preprocess does that)

# Preprocess dyadic covariates
# covariates start at t= 1 and end at t=7

all_ally <- preprocess(all_ally_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_ally)
sapply(all_ally, dim)

all_cont <- preprocess(all_cont_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_cont)
sapply(all_cont, dim)

all_trade <- preprocess(all_trade_nets, all_dep, 
                        lag = FALSE, covariate = TRUE)
str(all_trade)
sapply(all_trade, dim)


## Set vertex attributes

# get vertex names for the dependent networks
str(all_dep[[1]])
stateIDs1975 <- colnames(all_dep[[1]])
stateIDs1980 <- colnames(all_dep[[2]])
stateIDs1985 <- colnames(all_dep[[3]])
stateIDs1990 <- colnames(all_dep[[4]])
stateIDs1995 <- colnames(all_dep[[5]])
stateIDs2000 <- colnames(all_dep[[6]])
stateIDs2005 <- colnames(all_dep[[7]])

# turn dep into networks to be able to set vertex attributes
for (i in 1:length(all_dep)) {
  all_dep[[i]] <- network(all_dep[[i]])}
all_dep[[1]]

# For 1970 covariates (1975 dependent net)

pop1970 <- vcov1970$log_pop[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"pop", pop1970)

gdp1970 <- vcov1970$log_realgdp[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"gdp", gdp1970)

gdppc1970 <- vcov1970$log_rgdppc[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"gdppc", gdppc1970)

log_milex1970 <- vcov1970$log_milex[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"milex", log_milex1970)

log_milper1970 <- vcov1970$log_milper[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"milper", log_milper1970)

log_cinc1970 <- vcov1970$log_cinc[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"cinc", log_cinc1970)

nw1970 <- vcov1970$nw[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"nw", nw1970)

polity1970 <- vcov1970$polity[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"polity", polity1970)

dem_dum1970 <- vcov1970$dem_dum[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"dem_dum", dem_dum1970)

region1970 <- vcov1970$sub_region_n[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"region", region1970)

igo_hq1970 <- vcov1970$igo_hq[vcov1970$stateabb %in% stateIDs1975]
set.vertex.attribute(all_dep[[1]],"igo_hq", igo_hq1970)

# check
all_dep[[1]]


# For 1975 covariates (1980 dependent net)

pop1975 <- vcov1975$log_pop[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"pop", pop1975)

gdp1975 <- vcov1975$log_realgdp[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"gdp", gdp1975)

gdppc1975 <- vcov1975$log_rgdppc[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"gdppc", gdppc1975)

log_milex1975 <- vcov1975$log_milex[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"milex", log_milex1975)

log_milper1975 <- vcov1975$log_milper[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"milper", log_milper1975)

log_cinc1975 <- vcov1975$log_cinc[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"cinc", log_cinc1975)

nw1975 <- vcov1975$nw[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"nw", nw1975)

polity1975 <- vcov1975$polity[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"polity", polity1975)

dem_dum1975 <- vcov1975$dem_dum[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"dem_dum", dem_dum1975)

region1975 <- vcov1975$sub_region_n[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"region", region1975)

igo_hq1975 <- vcov1975$igo_hq[vcov1975$stateabb %in% stateIDs1980]
set.vertex.attribute(all_dep[[2]],"igo_hq", igo_hq1975)

# check
all_dep[[2]]



# For 1980 covariates (1985 dependent net)

pop1980 <- vcov1980$log_pop[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"pop", pop1980)

gdp1980 <- vcov1980$log_realgdp[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"gdp", gdp1980)

gdppc1980 <- vcov1980$log_rgdppc[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"gdppc", gdppc1980)

log_milex1980 <- vcov1980$log_milex[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"milex", log_milex1980)

log_milper1980 <- vcov1980$log_milper[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"milper", log_milper1980)

log_cinc1980 <- vcov1980$log_cinc[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"cinc", log_cinc1980)

nw1980 <- vcov1980$nw[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"nw", nw1980)

polity1980 <- vcov1980$polity[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"polity", polity1980)

dem_dum1980 <- vcov1980$dem_dum[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"dem_dum", dem_dum1980)

region1980 <- vcov1980$sub_region_n[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"region", region1980)

igo_hq1980 <- vcov1980$igo_hq[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[3]],"igo_hq", igo_hq1980)

# check
all_dep[[3]]


# For 1985 covariates (1990 dependent net)

pop1985 <- vcov1985$log_pop[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"pop", pop1985)

gdp1985 <- vcov1985$log_realgdp[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"gdp", gdp1985)

gdppc1985 <- vcov1985$log_rgdppc[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"gdppc", gdppc1985)

log_milex1985 <- vcov1985$log_milex[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"milex", log_milex1985)

log_milper1985 <- vcov1985$log_milper[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"milper", log_milper1985)

log_cinc1985 <- vcov1985$log_cinc[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"cinc", log_cinc1985)

nw1985 <- vcov1985$nw[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"nw", nw1985)

polity1985 <- vcov1985$polity[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"polity", polity1985)

dem_dum1985 <- vcov1985$dem_dum[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"dem_dum", dem_dum1985)

region1985 <- vcov1985$sub_region_n[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"region", region1985)

igo_hq1985 <- vcov1985$igo_hq[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[4]],"igo_hq", igo_hq1985)

# check
all_dep[[4]]


# For 1990 covariates (1995 dependent net)

pop1990 <- vcov1990$log_pop[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"pop", pop1990)

gdp1990 <- vcov1990$log_realgdp[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"gdp", gdp1990)

gdppc1990 <- vcov1990$log_rgdppc[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"gdppc", gdppc1990)

log_milex1990 <- vcov1990$log_milex[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"milex", log_milex1990)

log_milper1990 <- vcov1990$log_milper[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"milper", log_milper1990)

log_cinc1990 <- vcov1990$log_cinc[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"cinc", log_cinc1990)

nw1990 <- vcov1990$nw[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"nw", nw1990)

polity1990 <- vcov1990$polity[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"polity", polity1990)

dem_dum1990 <- vcov1990$dem_dum[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"dem_dum", dem_dum1990)

region1990 <- vcov1990$sub_region_n[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"region", region1990)

igo_hq1990 <- vcov1990$igo_hq[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[5]],"igo_hq", igo_hq1990)

# check
all_dep[[5]]


# For 1995 covariates (2000 dependent net)

pop1995 <- vcov1995$log_pop[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"pop", pop1995)

gdp1995 <- vcov1995$log_realgdp[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"gdp", gdp1995)

gdppc1995 <- vcov1995$log_rgdppc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"gdppc", gdppc1995)

log_milex1995 <- vcov1995$log_milex[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"milex", log_milex1995)

log_milper1995 <- vcov1995$log_milper[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"milper", log_milper1995)

log_cinc1995 <- vcov1995$log_cinc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"cinc", log_cinc1995)

nw1995 <- vcov1995$nw[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"nw", nw1995)

polity1995 <- vcov1995$polity[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"polity", polity1995)

dem_dum1995 <- vcov1995$dem_dum[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"dem_dum", dem_dum1995)

region1995 <- vcov1995$sub_region_n[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"region", region1995)

igo_hq1995 <- vcov1995$igo_hq[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[6]],"igo_hq", igo_hq1995)

# check
all_dep[[6]]


# For 2000 covariates (2005 dependent net)

pop2000 <- vcov2000$log_pop[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"pop", pop2000)

gdp2000 <- vcov2000$log_realgdp[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"gdp", gdp2000)

gdppc2000 <- vcov2000$log_rgdppc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"gdppc", gdppc2000)

log_milex2000 <- vcov2000$log_milex[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"milex", log_milex2000)

log_milper2000 <- vcov2000$log_milper[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"milper", log_milper2000)

log_cinc2000 <- vcov2000$log_cinc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"cinc", log_cinc2000)

nw2000 <- vcov2000$nw[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"nw", nw2000)

polity2000 <- vcov2000$polity[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"polity", polity2000)

dem_dum2000 <- vcov2000$dem_dum[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"dem_dum", dem_dum2000)

region2000 <- vcov2000$sub_region_n[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"region", region2000)

igo_hq2000 <- vcov2000$igo_hq[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[7]],"igo_hq", igo_hq2000)

# check
all_dep[[7]]


############################################ Estimate model #########################################################

set.seed(10)
# Model 4: democracy
model1a <- btergm(all_dep ~ edges + istar(2) + ostar(2) + mutual + triangles
                  +absdiff("polity")
                  +absdiff("gdppc")
                  +absdiff("milex")
                  +absdiff("nw")
                  +nodeicov("polity")
                  +nodeicov("gdppc")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)


################################################ Model2a ############################################################

### Load data

# Read in edgelists
el1980 <- read.csv("mod2_el1980.csv",stringsAsFactors=F)
el1985 <- read.csv("mod2_el1985.csv",stringsAsFactors=F)
el1990 <- read.csv("mod2_el1990.csv",stringsAsFactors=F)
el1995 <- read.csv("mod2_el1995.csv",stringsAsFactors=F)
el2000 <- read.csv("mod2_el2000.csv",stringsAsFactors=F)
el2005 <- read.csv("mod2_el2005.csv",stringsAsFactors=F)


# Read in dyadic covariates

# alliance networks
ally1980 <- read.csv("mod2_ally1980.csv",row.names=1)
ally1980 <- as.matrix(ally1980)
ally1985 <- read.csv("mod2_ally1985.csv",row.names=1)
ally1985 <- as.matrix(ally1985)
ally1990 <- read.csv("mod2_ally1990.csv",row.names=1)
ally1990 <- as.matrix(ally1990)
ally1995 <- read.csv("mod2_ally1995.csv",row.names=1)
ally1995 <- as.matrix(ally1995)
ally2000 <- read.csv("mod2_ally2000.csv",row.names=1)
ally2000 <- as.matrix(ally2000)
ally2005 <- read.csv("mod2_ally2003.csv",row.names=1)
ally2005 <- as.matrix(ally2005)

# contiguity matrix
cont1980 <- read.csv("mod2_cont1980.csv",row.names=1)
cont1980 <- as.matrix(cont1980)
cont1985 <- read.csv("mod2_cont1985.csv",row.names=1)
cont1985 <- as.matrix(cont1985)
cont1990 <- read.csv("mod2_cont1990.csv",row.names=1)
cont1990 <- as.matrix(cont1990)
cont1995 <- read.csv("mod2_cont1995.csv",row.names=1)
cont1995 <- as.matrix(cont1995)
cont2000 <- read.csv("mod2_cont2000.csv",row.names=1)
cont2000 <- as.matrix(cont2000)
cont2005 <- read.csv("mod2_cont2005.csv",row.names=1)
cont2005 <- as.matrix(cont2005)

# trade network
trade1980 <- read.csv("mod2_trade1980.csv",row.names=1)
trade1980 <- as.matrix(trade1980)
trade1985 <- read.csv("mod2_trade1985.csv",row.names=1)
trade1985 <- as.matrix(trade1985)
trade1990 <- read.csv("mod2_trade1990.csv",row.names=1)
trade1990 <- as.matrix(trade1990)
trade1995 <- read.csv("mod2_trade1995.csv",row.names=1)
trade1995 <- as.matrix(trade1995)
trade2000 <- read.csv("mod2_trade2000.csv",row.names=1)
trade2000 <- as.matrix(trade2000)
trade2005 <- read.csv("mod2_trade2005.csv",row.names=1)
trade2005 <- as.matrix(trade2005)

# Read in vertex covariates
vcov <- read.csv("mod2_vcov.csv",stringsAsFactors=F)
vcov1980 <- vcov[vcov$year==1980,]
vcov1985 <- vcov[vcov$year==1985,]
vcov1990 <- vcov[vcov$year==1990,]
vcov1995 <- vcov[vcov$year==1995,]
vcov2000 <- vcov[vcov$year==2000,]
vcov2005 <- vcov[vcov$year==2005,]


### Create network

# Initialize networks and assign labels

net1980 <- network.initialize(nrow(vcov1980))
network.vertex.names(net1980) <- vcov1980$stateabb

net1985 <- network.initialize(nrow(vcov1985))
network.vertex.names(net1985) <- vcov1985$stateabb

net1990 <- network.initialize(nrow(vcov1990))
network.vertex.names(net1990) <- vcov1990$stateabb

net1995 <- network.initialize(nrow(vcov1995))
network.vertex.names(net1995) <- vcov1995$stateabb

net2000 <- network.initialize(nrow(vcov2000))
network.vertex.names(net2000) <- vcov2000$stateabb

net2005 <- network.initialize(nrow(vcov2005))
network.vertex.names(net2005) <- vcov2005$stateabb


# Add in edges
net1980[as.matrix(el1980)] <- 1
net1985[as.matrix(el1985)] <- 1
net1990[as.matrix(el1990)] <- 1
net1995[as.matrix(el1995)] <- 1
net2000[as.matrix(el2000)] <- 1
net2005[as.matrix(el2005)] <- 1

# Define vertex attributes

# For 1980

set.vertex.attribute(net1980,"pop", vcov1980$log_pop)
set.vertex.attribute(net1980,"gdp", vcov1980$log_realgdp)
set.vertex.attribute(net1980,"gdppc", vcov1980$log_rgdppc)

set.vertex.attribute(net1980,"cinc", vcov1980$log_cinc)
set.vertex.attribute(net1980,"milex", vcov1980$log_milex)
set.vertex.attribute(net1980,"milper", vcov1980$log_milper)

set.vertex.attribute(net1980,"nw", vcov1980$nw)

set.vertex.attribute(net1980,"pts_inv", vcov1980$pts_inv)
set.vertex.attribute(net1980,"pts_dum", vcov1980$pts_dum)

set.vertex.attribute(net1980,"region", vcov1980$sub_region_n)
set.vertex.attribute(net1980,"igo_hq", vcov1980$igo_hq)


# For 1985

set.vertex.attribute(net1985,"pop", vcov1985$log_pop)
set.vertex.attribute(net1985,"gdp", vcov1985$log_realgdp)
set.vertex.attribute(net1985,"gdppc", vcov1985$log_rgdppc)

set.vertex.attribute(net1985,"cinc", vcov1985$log_cinc)
set.vertex.attribute(net1985,"milex", vcov1985$log_milex)
set.vertex.attribute(net1985,"milper", vcov1985$log_milper)

set.vertex.attribute(net1985,"nw", vcov1985$nw)

set.vertex.attribute(net1985,"pts_inv", vcov1985$pts_inv)
set.vertex.attribute(net1985,"pts_dum", vcov1985$pts_dum)

set.vertex.attribute(net1985,"region", vcov1985$sub_region_n)
set.vertex.attribute(net1985,"igo_hq", vcov1985$igo_hq)

# For 1990

set.vertex.attribute(net1990,"pop", vcov1990$log_pop)
set.vertex.attribute(net1990,"gdp", vcov1990$log_realgdp)
set.vertex.attribute(net1990,"gdppc", vcov1990$log_rgdppc)

set.vertex.attribute(net1990,"cinc", vcov1990$log_cinc)
set.vertex.attribute(net1990,"milex", vcov1990$log_milex)
set.vertex.attribute(net1990,"milper", vcov1990$log_milper)

set.vertex.attribute(net1990,"nw", vcov1990$nw)

set.vertex.attribute(net1990,"pts_inv", vcov1990$pts_inv)
set.vertex.attribute(net1990,"pts_dum", vcov1990$pts_dum)

set.vertex.attribute(net1990,"region", vcov1990$sub_region_n)
set.vertex.attribute(net1990,"igo_hq", vcov1990$igo_hq)

# For 1995

set.vertex.attribute(net1995,"pop", vcov1995$log_pop)
set.vertex.attribute(net1995,"gdp", vcov1995$log_realgdp)
set.vertex.attribute(net1995,"gdppc", vcov1995$log_rgdppc)

set.vertex.attribute(net1995,"cinc", vcov1995$log_cinc)
set.vertex.attribute(net1995,"milex", vcov1995$log_milex)
set.vertex.attribute(net1995,"milper", vcov1995$log_milper)

set.vertex.attribute(net1995,"nw", vcov1995$nw)

set.vertex.attribute(net1995,"pts_inv", vcov1995$pts_inv)
set.vertex.attribute(net1995,"pts_dum", vcov1995$pts_dum)

set.vertex.attribute(net1995,"region", vcov1995$sub_region_n)
set.vertex.attribute(net1995,"igo_hq", vcov1995$igo_hq)

# For 2000

set.vertex.attribute(net2000,"pop", vcov2000$log_pop)
set.vertex.attribute(net2000,"gdp", vcov2000$log_realgdp)
set.vertex.attribute(net2000,"gdppc", vcov2000$log_rgdppc)

set.vertex.attribute(net2000,"cinc", vcov2000$log_cinc)
set.vertex.attribute(net2000,"milex", vcov2000$log_milex)
set.vertex.attribute(net2000,"milper", vcov2000$log_milper)

set.vertex.attribute(net2000,"nw", vcov2000$nw)

set.vertex.attribute(net2000,"pts_inv", vcov2000$pts_inv)
set.vertex.attribute(net2000,"pts_dum", vcov2000$pts_dum)

set.vertex.attribute(net2000,"region", vcov2000$sub_region_n)
set.vertex.attribute(net2000,"igo_hq", vcov2000$igo_hq)

# For 2005

set.vertex.attribute(net2005,"pop", vcov2005$log_pop)
set.vertex.attribute(net2005,"gdp", vcov2005$log_realgdp)
set.vertex.attribute(net2005,"gdppc", vcov2005$log_rgdppc)

set.vertex.attribute(net2005,"cinc", vcov2005$log_cinc)
set.vertex.attribute(net2005,"milex", vcov2005$log_milex)
set.vertex.attribute(net2005,"milper", vcov2005$log_milper)

set.vertex.attribute(net2005,"nw", vcov2005$nw)

set.vertex.attribute(net2005,"pts_inv", vcov2005$pts_inv)
set.vertex.attribute(net2005,"pts_dum", vcov2005$pts_dum)

set.vertex.attribute(net2005,"region", vcov2005$sub_region_n)
set.vertex.attribute(net2005,"igo_hq", vcov2005$igo_hq)


# Check networks
net1980
net1985
net1990
net1995
net2000
net2005


### For the 1980-2005 period

# t=1 -> 1980
# t=2 -> 1985
# t=3 -> 1990
# t=4 -> 1995
# t=5 -> 2000
# t=6 -> 2005


# Combine dependent networks 

all_dip_nets <- list(as.matrix(net1980), as.matrix(net1985),
                     as.matrix(net1990), as.matrix(net1995), as.matrix(net2000), as.matrix(net2005))

# Combine dyadic covariates

all_ally_nets <- list(ally1980, ally1985, ally1990, ally1995, ally2000)
all_cont_nets <- list(cont1980, cont1985, cont1990, cont1995, cont2000)
all_trade_nets <- list(trade1980, trade1985, trade1990, trade1995, trade2000)
#all_ptrade_nets <- list(ptrade1980, ptrade1985, ptrade1990, ptrade1995, ptrade2000)


# check objects
class(all_dip_nets)
length(all_dip_nets)

# Use the preprocess function on the outcome network
# Tell it what missing values and structural zeros (instances in which tie impossible) look like

# estimation (dep nets) starts at t=2 and ends at t=6
all_dep <- preprocess(all_dip_nets,
                      lag = TRUE, covariate = FALSE, na = NA,
                      na.method = "fillmode", structzero = 10,
                      structzero.method = "remove") 

# Create lagged dependent network
all_mem <- preprocess(all_dip_nets, 
                      lag = TRUE, covariate = TRUE, memory = "autoregression",
                      na = NA, na.method = "fillmode", structzero = 10,
                      structzero.method = "remove")

# Check that everything looks right
length(all_dep)
sapply(all_dep, dim)
sapply(all_mem, dim)
# dim of the mem term need to be the same (preprocess does that)

# Preprocess dyadic covariates
# covariates start at t= 1 and end at t=5

all_ally <- preprocess(all_ally_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_ally)
sapply(all_ally, dim)

all_cont <- preprocess(all_cont_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_cont)
sapply(all_cont, dim)

all_trade <- preprocess(all_trade_nets, all_dep, 
                        lag = FALSE, covariate = TRUE)
str(all_trade)
sapply(all_trade, dim)


## Set vertex attributes

# get vertex names for the dependent networks
str(all_dep[[1]])
stateIDs1985 <- colnames(all_dep[[1]])
stateIDs1990 <- colnames(all_dep[[2]])
stateIDs1995 <- colnames(all_dep[[3]])
stateIDs2000 <- colnames(all_dep[[4]])
stateIDs2005 <- colnames(all_dep[[5]])

# turn dep into networks to be able to set vertex attributes
for (i in 1:length(all_dep)) {
  all_dep[[i]] <- network(all_dep[[i]])}
all_dep[[1]]


# For 1980 covariates (1985 dependent net)

pop1980 <- vcov1980$log_pop[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"pop", pop1980)

gdp1980 <- vcov1980$log_realgdp[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"gdp", gdp1980)

gdppc1980 <- vcov1980$log_rgdppc[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"gdppc", gdppc1980)

log_milex1980 <- vcov1980$log_milex[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"milex", log_milex1980)

log_milper1980 <- vcov1980$log_milper[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"milper", log_milper1980)

log_cinc1980 <- vcov1980$log_cinc[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"cinc", log_cinc1980)

nw1980 <- vcov1980$nw[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"nw", nw1980)

pts_inv1980 <- vcov1980$pts_inv[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"pts_inv", pts_inv1980)

pts_dum1980 <- vcov1980$pts_dum[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"pts_dum", pts_dum1980)

region1980 <- vcov1980$sub_region_n[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"region", region1980)

igo_hq1980 <- vcov1980$igo_hq[vcov1980$stateabb %in% stateIDs1985]
set.vertex.attribute(all_dep[[1]],"igo_hq", igo_hq1980)

# check
all_dep[[1]]


# For 1985 covariates (1990 dependent net)

pop1985 <- vcov1985$log_pop[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"pop", pop1985)

gdp1985 <- vcov1985$log_realgdp[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"gdp", gdp1985)

gdppc1985 <- vcov1985$log_rgdppc[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"gdppc", gdppc1985)

log_milex1985 <- vcov1985$log_milex[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"milex", log_milex1985)

log_milper1985 <- vcov1985$log_milper[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"milper", log_milper1985)

log_cinc1985 <- vcov1985$log_cinc[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"cinc", log_cinc1985)

nw1985 <- vcov1985$nw[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"nw", nw1985)

pts_inv1985 <- vcov1985$pts_inv[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"pts_inv", pts_inv1985)

pts_dum1985 <- vcov1985$pts_dum[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"pts_dum", pts_dum1985)

region1985 <- vcov1985$sub_region_n[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"region", region1985)

igo_hq1985 <- vcov1985$igo_hq[vcov1985$stateabb %in% stateIDs1990]
set.vertex.attribute(all_dep[[2]],"igo_hq", igo_hq1985)

# check
all_dep[[2]]


# For 1990 covariates (1995 dependent net)

pop1990 <- vcov1990$log_pop[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"pop", pop1990)

gdp1990 <- vcov1990$log_realgdp[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"gdp", gdp1990)

gdppc1990 <- vcov1990$log_rgdppc[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"gdppc", gdppc1990)

log_milex1990 <- vcov1990$log_milex[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"milex", log_milex1990)

log_milper1990 <- vcov1990$log_milper[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"milper", log_milper1990)

log_cinc1990 <- vcov1990$log_cinc[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"cinc", log_cinc1990)

nw1990 <- vcov1990$nw[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"nw", nw1990)

pts_inv1990 <- vcov1990$pts_inv[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"pts_inv", pts_inv1990)

pts_dum1990 <- vcov1990$pts_dum[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"pts_dum", pts_dum1990)

region1990 <- vcov1990$sub_region_n[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"region", region1990)

igo_hq1990 <- vcov1990$igo_hq[vcov1990$stateabb %in% stateIDs1995]
set.vertex.attribute(all_dep[[3]],"igo_hq", igo_hq1990)

# check
all_dep[[3]]


# For 1995 covariates (2000 dependent net)

pop1995 <- vcov1995$log_pop[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"pop", pop1995)

gdp1995 <- vcov1995$log_realgdp[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"gdp", gdp1995)

gdppc1995 <- vcov1995$log_rgdppc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"gdppc", gdppc1995)

log_milex1995 <- vcov1995$log_milex[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"milex", log_milex1995)

log_milper1995 <- vcov1995$log_milper[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"milper", log_milper1995)

log_cinc1995 <- vcov1995$log_cinc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"cinc", log_cinc1995)

nw1995 <- vcov1995$nw[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"nw", nw1995)

pts_inv1995 <- vcov1995$pts_inv[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"pts_inv", pts_inv1995)

pts_dum1995 <- vcov1995$pts_dum[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"pts_dum", pts_dum1995)

region1995 <- vcov1995$sub_region_n[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"region", region1995)

igo_hq1995 <- vcov1995$igo_hq[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[4]],"igo_hq", igo_hq1995)

# check
all_dep[[4]]


# For 2000 covariates (2005 dependent net)

pop2000 <- vcov2000$log_pop[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"pop", pop2000)

gdp2000 <- vcov2000$log_realgdp[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"gdp", gdp2000)

gdppc2000 <- vcov2000$log_rgdppc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"gdppc", gdppc2000)

log_milex2000 <- vcov2000$log_milex[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"milex", log_milex2000)

log_milper2000 <- vcov2000$log_milper[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"milper", log_milper2000)

log_cinc2000 <- vcov2000$log_cinc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"cinc", log_cinc2000)

nw2000 <- vcov2000$nw[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"nw", nw2000)

pts_inv2000 <- vcov2000$pts_inv[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"pts_inv", pts_inv2000)

pts_dum2000 <- vcov2000$pts_dum[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"pts_dum", pts_dum2000)

region2000 <- vcov2000$sub_region_n[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"region", region2000)

igo_hq2000 <- vcov2000$igo_hq[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[5]],"igo_hq", igo_hq2000)

# check
all_dep[[5]]


############################################ Estimate model #######################################################

set.seed(10)
# Model 5: human rights
model2a <- btergm(all_dep ~ edges + istar(2) + ostar(2) + mutual + triangles
                  +absdiff("pts_inv")
                  +absdiff("gdppc")
                  +absdiff("milex")
                  +absdiff("nw")
                  +nodeicov("pts_inv")
                  +nodeicov("gdppc")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)


################################################ Model 3a ##########################################################

### Load data

# Read in edgelists
el1995 <- read.csv("mod3_el1995.csv",stringsAsFactors=F)
el2000 <- read.csv("mod3_el2000.csv",stringsAsFactors=F)
el2005 <- read.csv("mod3_el2005.csv",stringsAsFactors=F)

# Read in dyadic covariates

# alliance networks
ally1995 <- read.csv("mod3_ally1995.csv",row.names=1)
ally1995 <- as.matrix(ally1995)
ally2000 <- read.csv("mod3_ally2000.csv",row.names=1)
ally2000 <- as.matrix(ally2000)
ally2005 <- read.csv("mod3_ally2003.csv",row.names=1)
ally2005 <- as.matrix(ally2005)

# contiguity matrix
cont1995 <- read.csv("mod3_cont1995.csv",row.names=1)
cont1995 <- as.matrix(cont1995)
cont2000 <- read.csv("mod3_cont2000.csv",row.names=1)
cont2000 <- as.matrix(cont2000)
cont2005 <- read.csv("mod3_cont2005.csv",row.names=1)
cont2005 <- as.matrix(cont2005)

# trade network
trade1995 <- read.csv("mod3_trade1995.csv",row.names=1)
trade1995 <- as.matrix(trade1995)
trade2000 <- read.csv("mod3_trade2000.csv",row.names=1)
trade2000 <- as.matrix(trade2000)
trade2005 <- read.csv("mod3_trade2005.csv",row.names=1)
trade2005 <- as.matrix(trade2005)

# Read in vertex covariates
vcov <- read.csv("mod3_vcov.csv",stringsAsFactors=F)
vcov1995 <- vcov[vcov$year==1995,]
vcov2000 <- vcov[vcov$year==2000,]
vcov2005 <- vcov[vcov$year==2005,]


### Create network

# Initialize networks and assign labels

net1995 <- network.initialize(nrow(vcov1995))
network.vertex.names(net1995) <- vcov1995$stateabb

net2000 <- network.initialize(nrow(vcov2000))
network.vertex.names(net2000) <- vcov2000$stateabb

net2005 <- network.initialize(nrow(vcov2005))
network.vertex.names(net2005) <- vcov2005$stateabb


# Add in edges
net1995[as.matrix(el1995)] <- 1
net2000[as.matrix(el2000)] <- 1
net2005[as.matrix(el2005)] <- 1

# Define vertex attributes

# For 1995

set.vertex.attribute(net1995,"pop", vcov1995$log_pop)
set.vertex.attribute(net1995,"gdp", vcov1995$log_realgdp)
set.vertex.attribute(net1995,"gdppc", vcov1995$log_rgdppc)

set.vertex.attribute(net1995,"cinc", vcov1995$log_cinc)
set.vertex.attribute(net1995,"milex", vcov1995$log_milex)
set.vertex.attribute(net1995,"milper", vcov1995$log_milper)

set.vertex.attribute(net1995,"nw", vcov1995$nw)

set.vertex.attribute(net1995,"hf_efiscore", vcov1995$hf_efiscore)
set.vertex.attribute(net1995,"ief_dum", vcov1995$ief_dum)

set.vertex.attribute(net1995,"region", vcov1995$sub_region_n)
set.vertex.attribute(net1995,"igo_hq", vcov1995$igo_hq)

# For 2000

set.vertex.attribute(net2000,"pop", vcov2000$log_pop)
set.vertex.attribute(net2000,"gdp", vcov2000$log_realgdp)
set.vertex.attribute(net2000,"gdppc", vcov2000$log_rgdppc)

set.vertex.attribute(net2000,"cinc", vcov2000$log_cinc)
set.vertex.attribute(net2000,"milex", vcov2000$log_milex)
set.vertex.attribute(net2000,"milper", vcov2000$log_milper)

set.vertex.attribute(net2000,"nw", vcov2000$nw)

set.vertex.attribute(net2000,"hf_efiscore", vcov2000$hf_efiscore)
set.vertex.attribute(net2000,"ief_dum", vcov2000$ief_dum)

set.vertex.attribute(net2000,"region", vcov2000$sub_region_n)
set.vertex.attribute(net2000,"igo_hq", vcov2000$igo_hq)

# For 2005

set.vertex.attribute(net2005,"pop", vcov2005$log_pop)
set.vertex.attribute(net2005,"gdp", vcov2005$log_realgdp)
set.vertex.attribute(net2005,"gdppc", vcov2005$log_rgdppc)

set.vertex.attribute(net2005,"cinc", vcov2005$log_cinc)
set.vertex.attribute(net2005,"milex", vcov2005$log_milex)
set.vertex.attribute(net2005,"milper", vcov2005$log_milper)

set.vertex.attribute(net2005,"nw", vcov2005$nw)

set.vertex.attribute(net2005,"hf_efiscore", vcov2005$hf_efiscore)
set.vertex.attribute(net2005,"ief_dum", vcov2005$ief_dum)

set.vertex.attribute(net2005,"region", vcov2005$sub_region_n)
set.vertex.attribute(net2005,"igo_hq", vcov2005$igo_hq)


# Check networks
net1995
net2000
net2005


### For the 1995-2005 period

# t=1 -> 1995
# t=2 -> 2000
# t=3 -> 2005


# Combine dependent networks 

all_dip_nets <- list(as.matrix(net1995), as.matrix(net2000), as.matrix(net2005))

# Combine dyadic covariates

all_ally_nets <- list(ally1995, ally2000)
all_cont_nets <- list(cont1995, cont2000)
all_trade_nets <- list(trade1995, trade2000)
#all_ptrade_nets <- list(ptrade1995, ptrade2000)


# check objects
class(all_dip_nets)
length(all_dip_nets)

# Use the preprocess function on the outcome network
# Tell it what missing values and structural zeros (instances in which tie impossible) look like

# estimation (dep nets) starts at t=2 and ends at t=3
all_dep <- preprocess(all_dip_nets,
                      lag = TRUE, covariate = FALSE, na = NA,
                      na.method = "fillmode", structzero = 10,
                      structzero.method = "remove") 

# Create lagged dependent network
all_mem <- preprocess(all_dip_nets, 
                      lag = TRUE, covariate = TRUE, memory = "autoregression",
                      na = NA, na.method = "fillmode", structzero = 10,
                      structzero.method = "remove")

# Check that everything looks right
length(all_dep)
sapply(all_dep, dim)
sapply(all_mem, dim)
# dim of the mem term need to be the same (preprocess does that)

# Preprocess dyadic covariates
# covariates start at t= 1 and end at t=5

all_ally <- preprocess(all_ally_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_ally)
sapply(all_ally, dim)

all_cont <- preprocess(all_cont_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_cont)
sapply(all_cont, dim)

all_trade <- preprocess(all_trade_nets, all_dep, 
                        lag = FALSE, covariate = TRUE)
str(all_trade)
sapply(all_trade, dim)


## Set vertex attributes

# get vertex names for the dependent networks
str(all_dep[[1]])
stateIDs2000 <- colnames(all_dep[[1]])
stateIDs2005 <- colnames(all_dep[[2]])

# turn dep into networks to be able to set vertex attributes
for (i in 1:length(all_dep)) {
  all_dep[[i]] <- network(all_dep[[i]])}
all_dep[[1]]


# For 1995 covariates (2000 dependent net)

pop1995 <- vcov1995$log_pop[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"pop", pop1995)

gdp1995 <- vcov1995$log_realgdp[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"gdp", gdp1995)

gdppc1995 <- vcov1995$log_rgdppc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"gdppc", gdppc1995)

log_milex1995 <- vcov1995$log_milex[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"milex", log_milex1995)

log_milper1995 <- vcov1995$log_milper[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"milper", log_milper1995)

log_cinc1995 <- vcov1995$log_cinc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"cinc", log_cinc1995)

nw1995 <- vcov1995$nw[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"nw", nw1995)

hf_efiscore1995 <- vcov1995$hf_efiscore[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"hf_efiscore", hf_efiscore1995)

ief_dum1995 <- vcov1995$ief_dum[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"ief_dum", ief_dum1995)

region1995 <- vcov1995$sub_region_n[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"region", region1995)

igo_hq1995 <- vcov1995$igo_hq[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"igo_hq", igo_hq1995)

# check
all_dep[[1]]


# For 2000 covariates (2005 dependent net)

pop2000 <- vcov2000$log_pop[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"pop", pop2000)

gdp2000 <- vcov2000$log_realgdp[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"gdp", gdp2000)

gdppc2000 <- vcov2000$log_rgdppc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"gdppc", gdppc2000)

log_milex2000 <- vcov2000$log_milex[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"milex", log_milex2000)

log_milper2000 <- vcov2000$log_milper[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"milper", log_milper2000)

log_cinc2000 <- vcov2000$log_cinc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"cinc", log_cinc2000)

nw2000 <- vcov2000$nw[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"nw", nw2000)

hf_efiscore2000 <- vcov2000$hf_efiscore[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"hf_efiscore", hf_efiscore2000)

ief_dum2000 <- vcov2000$ief_dum[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"ief_dum", ief_dum2000)

region2000 <- vcov2000$sub_region_n[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"region", region2000)

igo_hq2000 <- vcov2000$igo_hq[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"igo_hq", igo_hq2000)

# check
all_dep[[2]]


############################################ Estimate model ########################################################

set.seed(10)
# Model 6: economic freedom
model3a <- btergm(all_dep ~ edges + istar(2) + ostar(2) + mutual + triangles
                  +absdiff("hf_efiscore")
                  +absdiff("gdppc")
                  +absdiff("milex")
                  +absdiff("nw")
                  +nodeicov("hf_efiscore")
                  +nodeicov("gdppc")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)


texreg(c(model1a, model2a, model3a, model4e), digits=3, single.row=T, sideways=T)

